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We investigate the orbits of compact binary systems during the final inspiral period before co- 
alescence by integrating numerically the second-order post-Newtonian equations of motion. We 
include spin-orbit and spin-spin coupling terms, which, according to a recent study by Levin [J. 
Levin, Phys. Rev. Lett. 84, 3515 (2000)], may cause the orbits to become chaotic. To examine this 
claim, we study the divergence of initially nearby phase-space trajectories and attempt to measure 
the Lyapunov exponent 7. Even for systems with maximally spinning objects and large spin-orbit 
misalignment angles, we find no chaotic behavior. For all the systems we consider, we can place 
a strict lower limit on the divergence time iL = I/7 that is many times greater than the typical 
inspiral time, suggesting that chaos should not adversely affect the detection of inspiral events by 
T— ( ■ upcoming gravitational-wave detectors. 
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^ ■ For the current generation of ground-based gravitational-wave (GW) detectors, such as LIGO, one of the most 

^~5 promising sources is the final inspiral of two compact stars (black holes or neutron stars in a relativistic binary orbit) 

l/~j 0' PI- To detect such an event successfully, one must be able to match theoretical GW templates to experimental 

O^ \ data containing a great deal of instrumental noise. This "matched filtering" technique has the potential of greatly 
increasing the effective signal-to-noise ratio for the detector |l^, O] . Even to specify quantitative upper limits on inspiral 

'T'* ' events (a major objective of early LIGO runs), it is critical to have strong confidence in the physical accuracy of the 

J^ , templates Q. 

f-^ ' Recent work has suggested that the orbits of two rapidly spinning compact objects may be chaotic |q, o. The 

C^ ] presence of chaos on the inspiral timescale (typically ~ 100 s as the frequency sweeps up from ^^ 10 to 10 Hz) could 

f — significantly reduce the probability of detection, even with GW templates that include spin effects . In the matched 

^^ filtering technique, if the signal and template are off by as little as half a wavelength (over ~ 10 cycles), the event 

^"^ I could be missed because of destructive interference B. In order to quantify this threat, we attempt to measure 

^^ ' Lyapunov divergence times for a broad sample of initial conditions. 

O As in the Newtonian two-body problem, a relativistic binary system can be expressed in terms of a reduced mass 

O"^ /i = mim2/{mi + 7712) orbiting a fixed mass m — nii + 7712 with a separation r = r 1 — r2. We adopt units with 

^ G = c = m = 1 and use the same notation and post-Newtonian (PN) equations of motion as described by Kidder 

f^ [g[, expanded to second order (2PN) terms in {v/cf' and Gm/r, and including spin-orbit (SO) and spin-spin (SS) 

►;>. \ effects as well as the 2.5PN radiation reaction (RR) term. The individual spins precess because of frame-dragging 

and the Lens-Thirring effect according to Si =(\i x S^. Here 7 = 1,2, Si is the spin of the compact object, and ^i is 
5^ ' a variable axis of precession as defined in eqn. (2.4) of M. The magnitude of Cli is the instantaneous spin precession 
Cd ' frequency. We also define a mass ratio /3 = 7712/7771, and a spin-orbit misalignment angle 9i for each object as the 
angle between the spin vector Si and the Newtonian angular momentum L^r = /^(r x ?). The relative position vector 
f evolves according to a second-order ordinary differential equation of the type r = apN + aso + ass + aRR. The 
full expressions for these terms can be found in Refs. 1^ and 0. 

We integrate these equations numerically in double precision using a 5th-order Runge-Kutta scheme with an adaptive 
time step. The robust nature of the Runge-Kutta algorithm makes it particularly attractive for measuring exponential 
divergence of nearby trajectories. Indeed, since the Runge-Kutta integration can only introduce errors that grow at 
a polynomial rate, any exponential divergence, if it occurs, will rapidly dominate the evolution and should be easily 
distinguishable from numerical effects |9[|. To quantify chaotic behavior in a dynamical system, the equations of 
motion must be conservative in the sense that there should be no dissipative terms that could act as attractors in 
phase space (which would eliminate the possibility of the system being formally chaotic; see |10| ). For this reason, 
when integrating the equations of motion to calculate Lyapunov exponents, the radiation reaction terms are not 
included. This also allows for the option of very long integrations to test whether the system is formally chaotic (but 
on a time scale much longer than the actual inspiral time). 
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Historically, an important tool for identifying chaos in celestial mechanics has been the Poincare surface of section 
|_l|. By plotting the position in phase space only at certain values of independent coordinate variables, one can 
reduce a complicated four-dimensional (4-D) phase space trajectory to a 2-D scatter plot. Conserved quantities such 
as energy or angular momentum generally are held constant for different initial conditions plotted on a single section. 
For many systems, some initial conditions behave regularly, producing 1-D curves as if there were additional integrals 
of the motion, while others fill out 2-D regions of the section. This spreading away from the invariant curves is 
evidence for chaotic behavior in the system. However, this method is only successful in reducing the phase space 
by one additional dimension, so for systems with [# of degrees of freedom] — [# of integrals of motion] greater than 
two, the surface of section technique is not very useful for identifying chaos |jl3]. For higher-dimensional systems, the 
projection of non-chaotic orbits onto a two-dimensional section in general will not generate confined curves, and thus 
a spreading of points in this pseudo-section is not necessarily an indication of chaos. For the problem of two spinning 
compact objects, the phase space is 10-D (3 from r, 3 from r, and 2 each from the spin vectors Si and S2, which can 
precess in arbitrary orientations but maintain constant magnitude [p3|) while there are only 4 obvious constraints, 
corresponding to angular momentum and energy conservation due to the invariance of the PN Lagrangian under 
rotations and time translations S. This almost guarantees that the projection of the trajectory onto a 2-D section 
of phase space will not be constrained to a 1-D curve, so regular behavior could easily be misinterpreted as chaos (cf. 

il). 

A more quantitative method for identifying and measuring chaos is to calculate the maximum Lyapunov exponent, 

defined as the divergence rate between initially nearby trajectories: 

, , 1 , f dX{t)\ 
-f(t) = - n H- ■ (1) 

Here the difference dX between two points in phase space is simply the Cartesian distance between the dimensionless 
12-component coordinate vectors [f, r. Si, S2] and [f, r', S'l, S2] of two nearby trajectories. In our numerical integra- 
tions, the initial separation is a small displacement in phase space with a random orientation and a magnitude of 
dX{0) — 10^^°. The two trajectories are then integrated forward in time, recording the separation dX{t), from which 
7 is computed. In chaotic systems, the divergence will be exponential in time with a roughly constant (positive) 
exponent: dX(t) ^ dX (0) exjpi^t) . To be precise, there are actually many different Lyapunov exponents, one for 
each dimension of phase space. The maximum Lyapunov exponent is automatically selected because, like a classical 
eigenvector problem, any vector (such as the random initial displacement vector) when multiplied repeatedly by the 
same matrix will grow fastest in the direction corresponding to the largest eigenvalue. Similarly, for a chaotic system, 
any random initial displacement in phase space is expected to grow as fast as the maximum Lyapunov exponent. 
However, for a regular system, the divergence generally grows linearly or at most as a power law in time. In that 
case, the Lyapunov exponent 7 approaches zero for large t. We define the Lyapunov time ^l = I/7 as the time scale 
on which nearby trajectories separate by a factor of e (the "e-folding time"). While it is difficult, on the basis of 
numerical integrations, to claim that a system is categorically non-chaotic on all time scales, we can make the more 
precise claim that on a particular time scale the system shows no chaotic behavior, that is, we can set a lower limit on 
tL. For the problem of coalescing compact binaries, the relevant time scale is the inspiral time tjnsp, so if tjnsp ■C tL, 
chaos will not affect the dynamics. 

To establish the validity of our numerical approach, we have re-examined the problem studied by Suzuki and 
Maeda [[l4| , where both chaotic and regular trajectories have clearly been identified for super-maximally spinning test 
particles (^2 » ml' w ^ /x) orbiting around a Schwarzschild black hole {Si = 0). While physically unreahstic, these 
conditions introduce no mathematical singularities into the equations of motion and are therefore fine examples for 
identifying chaotic or regular behavior in this dynamical system. Here we show numerical results for two different sets 
of initial conditions: both begin with r/m = 4.5, /3 — 10~^, S'l = 0, ^2 = 1.4 x 10^ rn^, and spin-orbit misalignment 
angle O2 = 7r/4, but one trajectory has an initial Newtonian orbital angular momentum of Ljq /m — 1.485 /^ while 
the other begins with Ljq jm = 1.53 /x. In Fig. 1 we show the measured value of 7 as a function of time (measured in 
orbital periods). The orbit with Ljq /m = 1.485 ^ (a) is clearly chaotic with 7 ~ 0.1 (corresponding to a divergence 
time scale of ^l ~ 10 orbital periods). The other trajectory (b) exhibits very different qualitative behavior, with no 
evidence for chaos even on the time scale of thousands of orbits. Our results for these and many other trajectories 
that we have computed agree well with Suzuki and Maeda's results, confirming that both chaotic and regular orbits 
exist for rapidly spinning test particles around a Schwarzschild black hole. 

Another method for distinguishing between chaotic and quasi-periodic orbits is to look for certain qualitative 
features in the power spectrum of one of the dynamical variables ]9|, 15|. For the compact binary problem, the obvious 
choice is to analyze the spectrum of the GW signal, defined as the squared amplitude of the Fourier transform of 
the GW strain h{t). With RR turned off, the spectrum should exhibit sharp lines for quasi-periodic orbits, but 
broad-band noise for chaotic orbits. Including only the leading quadrupole radiation terms, the components of the 



transverse-traceless GW tensor observed from the z direction are 

4/i / 2 "m 



h+ = 



D 



{-I ~ ^-') - (2) 



^x = -^ [v^Vy - -^xyj , (3) 

where D is the distance to the source and Vi — Xi [R 16| . Figure 2 shows the power spectra in h^ (t) for the same two 
test-particle orbits as in Fig. 1. Both spectra show peaks around /gw — 25 Hz, the fundamental quadrupolc frequency 
(twice the orbital frequency) for these initial conditions. As expected, the chaotic orbit (a) produces primarily broad- 
band noise with few distinguishable features, whereas the spectrum for the quasi-periodic orbit (b) shows many sharp 
lines. 

Having established the ability of our method to distinguish between regular and chaotic trajectories described by 
the PN equations of motion, we now apply the same techniques to the astrophysically relevant systems expected to be 
detectable by ground-based laser interferometers. Here we will show results for an illustrative case of two maximally 
spinning IOMq black holes {Si = mf) in an initially circular orbit with spin-orbit misalignment angles of 9i ~ 38° 
and 02 = 70°. To measure the Lyapunov exponent at different points during the final inspiral, 4 different trajectories 
were integrated, corresponding to GW frequencies of 10, 40, 100, and 400 Hz, emitted at orbital separations r/m ~ 
50, 20, 10, and 5, respectively |0|. For each separation, we then integrate the PN equations of motion (with RR 
turned off) over a period much longer than the inspiral time. If no exponential divergence is observed, we can safely 
conclude that the system is, for all practical purposes, not chaotic. 

The measured values for 7(i) are plotted in Fig. 3 for each of the selected stages of the inspiral. For any power-law 
divergence in phase space with dX{t) ^ t", on a log- log plot of 7 versus time the slope is dln7/dlnt = 1/ Int — 1, so 
that for large times the curve should be linear with slope -1 [as we see in Fig. 1 (b)]. All plots in Fig. 3 clearly show 
this behavior, characteristic of regular, non-chaotic orbits. Also shown for comparison in Fig. 3 are the inspiral times 
corresponding to each separation. In all four cases, the lower limit on ^l far exceeds tinsp- 

Just as in the test particle example presented above, we also calculate the gravitational wave power spectrum for 
this black hole binary. Here the results are even more striking: Fig. 4 shows the h^{t) power spectrum at the point 
in the inspiral where the GW frequency is 100 Hz. The spectrum has a few very sharp lines with no evidence for 
broad-band noise. The strongest features correspond to the fundamental GW quadrupole and spin-orbit precessional 
frequencies (100 Hz and 7 Hz, respectively) and their harmonics. The spectra from the other points in the inspiral 
show the same qualitative features, all indicative of quasi-periodic motion. 

In addition to this prototypical binary black hole system, we have investigated the behavior of a large number of 
other systems at different stages throughout the inspiral. Varying the binary mass ratio, spin magnitudes (always with 
Si < mf), misalignment angles, eccentricity, and initial separation over wide ranges, we consistently find the same 
regular, non-chaotic behavior for all trajectories. In particular, we have calculated orbits for the same high-eccentricity 
systems considered by Levin in Ref. H and again measure only linear divergence of nearby initial conditions, finding 
iinsp ^ ^L iQ- We have also looked at many binary systems in which either one or both objects are solar-mass 
neutron stars. Because of their smaller mass (and thus slower radiation loss), these systems can spend significantly 
longer times in the inspiral band, completing tens of thousands of orbits before merger. Yet, even on the longest 
relevant time scales (iinsp ^ 10^ s), we again find no chaotic behavior. 

Recent work in binary star evolution [|l9[ suggests that systems with large spin misalignment angles are likely to 
be astrophysically relevant sources, giving rise to complicated nonlinear SO and SS interactions. It is possible that 
there are undiscovered narrow regions in phase space where specific resonant conditions could give rise to chaotic 
behavior, yet even if a binary system passes through them during inspiral, the time spent in these resonant bands 
will most likely be much shorter than the Lyapunov time. Even though the Lyapunov exponent approaches zero 
for all the systems we have considered, their orbits can still appear quite irregular, especially when both objects are 
spinning rapidly pj. As mentioned above, this can lead to the spreading of the surface of section away from a 1-D 
curve, but in dynamical systems with many degrees of freedom, this spreading alone does not imply the presence 
of chaos. Nevertheless, the detection of GW signals from binaries with rapidly spinning components may be very 
challenging. Systems with similar initial conditions still may produce waveforms that look quite different even without 
exponential divergence. Orbital precession will add many complicated features to the basic inspiral waveforms. The 
current template database will most likely need to be extended (perhaps by orders of magnitude) to include at least 
a portion of parameter space where the signals are significantly modified by spin effects poj. Furthermore, for the 
late stages of inspiral, the 2.5PN expansion becomes invalid and one should include higher order PN terms (or even 
perform calculations in full general relativity pl|). While these terms introduce additional features in the dynamics, 
at this late stage the orbit is certainly decaying too rapidly for the final GW signal to be affected by formal chaotic 
behavior. 
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FIG. 1: Lyapunov exponent 7 for test-particle orbits around a black hole. The orbits in case (a) are chaotic, with 7 approaching 
a positive value corresponding to a divergence time of t^ ~ 10 orbital periods. The regular orbits in case (b) diverge linearly 
in time so that 7(f) -^ 0. 
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(b) Quasi — Period Orbit 
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FIG. 2: Power spectra of the GW signal h+{t) for the same two cases as in Fig. 1. The chaotic system in (a) produces broadband 
noise while the quasi-periodic orbit (b) exhibits sharp spectral lines. The GW frequency is given in units of m^g Hz, where 
niu) = m/lOMQ. 
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FIG. 3: Lyapunov exponent 7(i) for different stages in the inspiral of two maximally spinning WMq black holes: /gw ~ 10 Hz, 
r/m = 47.25, tinsp = 69 s; /gw = 40 Hz, r/m = 18.75, iinsp = 7.5 s; /gw = 100 Hz, r/m = 9.2, iin^p = 0.2 s; /gw = 400 Hz, 
r/m = 4.0, iinsp ~ 0.01 s. No evidence for chaos is seen, with II ~ I/7 S> tinsp in all cases. 
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FIG. 4: Power spectrum of the GW signal h+{t) calculated from the initial conditions of Fig. 3 at the stage in the inspiral 

Gw = 100 Hz. Tl 
that the orbit is regular. 



where /gw = 100 Hz. The frequency is in units off m2(/ Hz, where m2o ~ rn/20MQ. The sharp lines in the spectrum confirm 



